



full.data.models<-function(data=NULL, dv=NULL, int=F,  short.data=F, short.data.cut=NULL,  model.data=NULL,  vc=NULL, gen.vc=F, std=F){
	
	

	
	
if(int==F){	
	##generate lagged DV

		
#CONSTRUCT LAGGED DV
		
form<-paste(dv,"~","dayno")
form2<-as.formula(form)		
means<-summaryBy(form2, data=data, FUN=mean, na.rm=T)##get daily hit rates
means$lag.hit.new<-NA
means<-means[order(means$dayno, decreasing=F ),]
head(means)

for(i in 2:length(means$dayno)){
	
	means$lag.hit.new[i]<-means[i-1,2]
	
}


head(means)



data<-merge(data, means, by="dayno")

		
if(short.data==T){

	data<-subset(data, data$dayno>=(memo-short.data.cut) & data$dayno<(memo+short.data.cut))
	
}
	
	
se1<-rep(NA, 8)
names.models<-c("diff.means","diff.means+controls","linear","linear+controls","squared","squared+controls","cubic","cubic+controls")
names(se1)<-names.models

##diff in means
m1<-lm(data[,dv]~period, data=data)
se1[1]<-summary(m1)$coefficients["periodpost-memo","Std. Error"]
m1a<-lm(data[,dv]~period+factor(year)+factor(month)+factor(weekday)+lag.hit.new, data=data)
se1[2]<-summary(m1a)$coefficients["periodpost-memo","Std. Error"]

##local linear
m2<-lm(data[,dv]~period+distance+period*distance, data=data)
se1[3]<-summary(m2)$coefficients["periodpost-memo","Std. Error"]
m2a<-lm(data[,dv]~period+distance+period*distance+factor(year)+factor(month)+factor(weekday)+lag.hit.new, data=data)
se1[4]<-summary(m2a)$coefficients["periodpost-memo","Std. Error"]

##second order poly
m3<-lm(data[,dv]~period+distance+period*distance+I(distance^2)+period*I(distance^2), data=data)
se1[5]<-summary(m3)$coefficients["periodpost-memo","Std. Error"]
m3a<-lm(data[,dv]~period+distance+period*distance+I(distance^2)+period*I(distance^2)+factor(year)+factor(month)+factor(weekday)+lag.hit.new, data=data)
se1[6]<-summary(m3a)$coefficients["periodpost-memo","Std. Error"]


##cubic
m4<-lm(data[,dv]~period+distance+period*distance+I(distance^2)+period*I(distance^2)+I(distance^3)+period*I(distance^3), data=data)
se1[7]<-summary(m4)$coefficients["periodpost-memo","Std. Error"]

m4a<-lm(data[,dv]~period+distance+period*distance+I(distance^2)+period*I(distance^2)+I(distance^3)+period*I(distance^3)+factor(year)+factor(month)+factor(weekday)+lag.hit.new, data=data)
se1[8]<-summary(m4a)$coefficients["periodpost-memo","Std. Error"]

models2<-list(m1, m1a, m2, m2a, m3, m3a, m4, m4a)

print("models done")

##load HAC SEs
coefs<-NA
ses<-NA
N<-NA
se.final<-NA



if(!is.null(vc)){
load(vc)
se2<-NA
for(j in 1:length(vc)){
se2[j]<-sqrt(vc[[j]]["periodpost-memo","periodpost-memo"])
}
}

if(std==T){
se2<-se1
}


if(gen.vc==T){
vc<-list(NA)

for(k in 1:length(models2)){
	vc[[k]]<-vcovHAC(models2[[k]])	
}
se2<-NA
for(j in 1:length(vc)){
se2[j]<-sqrt(vc[[j]]["periodpost-memo","periodpost-memo"])
}
}



for(i in 1:length(models2)){
	if(se1[i]>se2[i]){se.final[i]<-se1[i]}
	if(se1[i]<se2[i]){se.final[i]<-se2[i]}
	if(se1[i]==se2[i]){se.final[i]<-se1[i]}

}

for(i in 1:length(models2)){	
	coefs[i]<-models2[[i]]$coefficients["periodpost-memo"]
	ses[i]<-se.final[i]
	N[i]<-nrow(models2[[i]]$model)
}

sig<-NA
for(i in 1:length(se.final)){
	
	lb<-coefs[i]-1.96*se.final[i]
	ub<-coefs[i]+1.96*se.final[i]
	if(lb*ub>0){sig[i]<-1}
	if(lb*ub<0){sig[i]<-0}
	
	
}


r<-rbind.data.frame(coefs, ses,N)
r[1:2,]<-round(r[1:2,], digits=3)
for(i in 1:ncol(r)){
	
	if(sig[i]==1){r[2,i]<-paste("(",r[2,i],")*",sep="")}
	if(sig[i]==0){r[2,i]<-paste("(",r[2,i],")",sep="")}
	
}
library(xtable)
rownames(r)<-c("tau","","N")
colnames(r)<-c("dm","dmc","l","lc","s","sc","cubic","cubicc")
table<-xtable(r, digits=2,include.rownames=FALSE)

results<-list(coefs, se.final, r, table)
names(results)<-c("coefs","se.final","r","table")
return(results)


}


if(int==T){
	load(model.data)##load conventional SEs
	load(vc)##load robust HAC SEs
	co<-NA
	coefs<-list(co, co, co, co, co, co, co,co)
	N<-NA
	se1<-NA
	se2<-NA
	sef<-NA
	se.final<-list(sef, sef, sef, sef, sef, sef, sef, sef)
	
	
	for(j in 1:length(models2)){
	coefs[[j]][1]<-models2[[j]]$coefficients["periodpost-memo"]
	coefs[[j]][2]<-models2[[j]]$coefficients["periodpost-memo"]+models2[[j]]$coefficients["periodpost-memo:int"]
	coefs[[j]][3]<-models2[[j]]$coefficients["periodpost-memo:int"]

    vcs<-vcov(models2[[j]])

	se1[1]<-summary(models2[[j]])$coefficients["periodpost-memo","Std. Error"]
	se1[2]<-sqrt(vcs["periodpost-memo","periodpost-memo"] + vcs["periodpost-memo:int","periodpost-memo:int"] + 2*vcs["periodpost-memo","periodpost-memo:int"])
	se1[3]<-summary(models2[[j]])$coefficients["periodpost-memo:int","Std. Error"]

	se2[1]<-sqrt(vc[[j]]["periodpost-memo","periodpost-memo"])	
	se2[2]<-sqrt(vc[[j]]["periodpost-memo","periodpost-memo"] + vc[[j]]["periodpost-memo:int","periodpost-memo:int"] + 2*vc[[j]]["periodpost-memo","periodpost-memo:int"])
	se2[3]<-sqrt(vc[[j]]["periodpost-memo:int","periodpost-memo:int"])
	
	
	
	for(f in 1:3){
	if(se1[f]>se2[f]){se.final[[j]][f]<-se1[f]}
	if(se1[f]<se2[f]){se.final[[j]][f]<-se2[f]}
	N[j]<-nrow(models2[[j]]$model)
	}
	}


sig1<-NA
sig<-list(sig1, sig1, sig1, sig1, sig1, sig1, sig1, sig1)
for(j in 1:length(se.final)){
	
	for(i in 1:length(se.final[[j]])){
	lb<-coefs[[j]][i]-1.96*se.final[[j]][i]
	ub<-coefs[[j]][i]+1.96*se.final[[j]][i]
	if(lb*ub>0){sig[[j]][i]<-1}
	if(lb*ub<0){sig[[j]][i]<-0}
	
	}
}

r<-cbind.data.frame(rep(NA, 7))
for(j in 1:length(se.final)){
temp<-rbind.data.frame(coefs[[j]][1], se.final[[j]][1],coefs[[j]][2], se.final[[j]][2],coefs[[j]][3], se.final[[j]][3],N[j])
r<-cbind.data.frame(r, temp)
}
r<-r[,-1]
colnames(r)<-1:8
r<-round(r, digits=3)
r

for(i in 1:ncol(r)){
	if(sig[[i]][1]==1){r[2,i]<-paste("(",r[2,i],")*",sep="")}
	if(sig[[i]][1]==0){r[2,i]<-paste("(",r[2,i],")",sep="")}
	if(sig[[i]][2]==1){r[4,i]<-paste("(",r[4,i],")*",sep="")}
	if(sig[[i]][2]==0){r[4,i]<-paste("(",r[4,i],")",sep="")}
	if(sig[[i]][3]==1){r[6,i]<-paste("(",r[6,i],")*",sep="")}
	if(sig[[i]][3]==0){r[6,i]<-paste("(",r[6,i],")",sep="")}
}
r

library(xtable)
rownames(r)<-c("no.group","se1","group","se2","diff","se3","N")
colnames(r)<-c("dm","dmc","l","lc","s","sc","cubic","cubicc")
table<-xtable(r, digits=2,include.rownames=FALSE)

results<-list(coefs, se.final, r, table)
names(results)<-c("coefs","se.final","r","table")
return(results)
}

}

	
	
	
	
	
	
	
	
	
